Dealing with many zeros
Count data are often modeled with Poisson or Negative Binomial models.
These assume zeros occur naturally from the count process.
But in real data, we often see too many zeros → zero inflation.
Examples:
Both handle excess zeros but use different assumptions about how zeros are generated.
Zero-Inflated Models
Hurdle Models
(Occurrence Models)
Both handle excess zeros but use different assumptions about how zeros are generated.
Zero-Inflated Models
Hurdle Models
Both handle excess zeros but use different assumptions about how zeros are generated:
Zero-Inflated Models
Hurdle Models
Assume two processes:
A binary process → decides if the observation is a certain zero.
A count process → generates counts, including zeros.
\[ \begin{aligned} \text{P}(Y=0) & = \pi + (1-\pi)f(0) & \\ \text{P}(Y=y) & = (1-\pi)f(y)& y=1,2,\dots \end{aligned} \]
Note In zero-inflated models \(f(y)\) is typically discrete (ex Poisson, Binomial, …)
Also have two parts, but different logic:
A binary model determines if the observation crosses the hurdle (i.e. zero vs. positive).
A truncated count model models only positive counts.
Zeros come only from the first process.
\[ \begin{aligned} \text{P}(Y=0) & = \pi & \\ \text{P}(Y=y) & = \frac{f(y)}{1-f(0)}& y=1,2,\dots \end{aligned} \]
Note In hurdle models \(f_Y(y)\) can also be continuous (f.ex Gamma). In these cases: \[ \begin{aligned} \text{P}(Y=0) & = \pi & \\ f_Y(y) & = f(y)& y>0 \end{aligned} \] where \(f(y)\) is scaled so that it integrates to \(1-\pi\)
| Feature | Zero-Inflated | Hurdle |
|---|---|---|
| Zeros come from | Two sources (inflation + count) | One source (hurdle only) |
| Positive part | Includes zeros | Truncated at zero |
| Can use continuous distributions? | ❌ Typically count only | ✅ Yes (Gamma, Lognormal) |
inlabru (Type 1)zeroinflatedpoisson1)zeroinflatedbinomial1)zeroinflatednbinomial1)zeroinflatedbinomial1)To get details about the distributions you can type
The probability of the inflation process \(\pi\) is a hyperparameter therefore cannot be modeled with covariates
Number of fishes caught by fishermen at a state park.
nofish persons child count
1 1 1 0 0
2 0 1 0 0
3 0 1 0 0
\[ \begin{eqnarray} \text{Model 1: }\ & Y\sim \text{Poisson}(\lambda)\\ \text{Model 2: }\ & Y\sim \text{NegBinomial}(n,p)\\ \text{Linear predictor: }\ & \eta = \beta_0 + \beta_1x_1 + \beta_2x_2 \end{eqnarray} \]
\[ \begin{eqnarray} \text{Model 1: }\ & Y\sim \text{Poisson}(\lambda)\\ \text{Model 2: }\ & Y\sim \text{NegBinomial}(n,p)\\ \text{Linear predictor: }\ & \eta = \beta_0 + \beta_1x_1 + \beta_2x_2 \end{eqnarray} \]
Negative Binomial distribution
\[ \text{Prob}(Y=y) = \frac{\Gamma(y+n)}{\Gamma(n)\Gamma(y+1)}p^n(1-p)^y \] with \[ E(Y) = \mu = n\frac{1-p}{p} \qquad \text{Var}(Y) = \mu(1+\frac{\mu}{n}) \] and linear predictor \(\eta = log(\mu)\)
The Model \[ \begin{aligned} \text{Model 1 }: y_i|\eta_i & \sim \pi1_{(y=0)} + (1-\pi)\text{Pois}(\lambda(\eta_i))\\ \text{Model 2 }: y_i|\eta_i & \sim \pi1_{(y=0)} + (1-\pi)\text{NegBin}(n(\eta_i),p(\eta_i))\\ \text{Linear Predictor }: \eta_i & = \color{red}{\boxed{\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i}}} \end{aligned} \]
The code
# define model components
cmp = ~ Intercept(1) + fishing(nofish, model = "linear") +
persons(persons, model = "linear") +
child(child, model = "linear")
# define model predictor
formula = count ~ .
# build the observation model
# Poisson model
lik_pois = bru_obs(formula,
family = "zeroinflatedpoisson1",
data = zinb)
# Negative binomial model
lik_nbin = bru_obs(formula,
family = "zeroinflatednbinomial1",
data = zinb)
# fit the model
fit_pois = bru(cmp, lik_pois)
fit_nbin = bru(cmp, lik_nbin)The Model \[ \begin{aligned} \text{Model 1 }: \color{red}{\boxed{y_i|\eta_i}} & \color{red}{\boxed{\sim \pi1_{(y=0)} + (1-\pi)\text{Pois}(\lambda(\eta_i))}}\\ \text{Model 2 }: y_i|\eta_i & \sim \pi1_{(y=0)} + (1-\pi)\text{NegBin}(n(\eta_i),p(\eta_i))\\ \text{Linear Predictor }: \eta_i & = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} \end{aligned} \]
The code
# define model components
cmp = ~ Intercept(1) + fishing(nofish, model = "linear") +
persons(persons, model = "linear") +
child(child, model = "linear")
# define model predictor
formula = count ~ .
# build the observation model
# Poisson model
lik_pois = bru_obs(formula,
family = "zeroinflatedpoisson1",
data = zinb)
# Negative binomial model
lik_nbin = bru_obs(formula,
family = "zeroinflatednbinomial1",
data = zinb)
# fit the model
fit_pois = bru(cmp, lik_pois)
fit_nbin = bru(cmp, lik_nbin)The Model \[ \begin{aligned} \text{Model 1 }: y_i|\eta_i & \sim \pi1_{(y=0)} + (1-\pi)\text{Pois}(\lambda(\eta_i))\\ \text{Model 2 }: \color{red}{\boxed{y_i|\eta_i}} & \color{red}{\boxed{\sim \pi1_{(y=0)} + (1-\pi)\text{NegBin}(n(\eta_i),p(\eta_i))}}\\ \text{Linear Predictor }: \eta_i & = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} \end{aligned} \]
The code
bru_options_set(control.compute = list(dic = T, waic = T, mlik = T))
# define model components
cmp = ~ Intercept(1) + fishing(nofish, model = "linear") +
persons(persons, model = "linear") +
child(child, model = "linear")
# define model predictor
formula = count ~ .
# build the observation model
# Poisson model
lik_pois = bru_obs(formula,
family = "zeroinflatedpoisson1",
data = zinb)
# Negative binomial model
lik_nbin = bru_obs(formula,
family = "zeroinflatednbinomial1",
data = zinb)
# fit the model
fit_pois = bru(cmp, lik_pois)
fit_nbin = bru(cmp, lik_nbin)[1] "Poisson Model"
mean sd 0.025quant 0.975quant
Intercept 0.22 0.15 -0.08 0.51
fishing -0.90 0.12 -1.13 -0.67
persons 0.65 0.05 0.56 0.74
child -0.95 0.10 -1.15 -0.76
[1] "Negative Binomial Model"
mean sd 0.025quant 0.975quant
Intercept -0.80 0.31 -1.42 -0.18
fishing -0.59 0.25 -1.08 -0.10
persons 0.94 0.12 0.71 1.17
child -1.65 0.19 -2.03 -1.27
Hyperparameters
[1] "Poisson Model"
mean sd
zero-probability parameter for zero-inflated poisson_1 0.48 0.04
[1] "Negative Binomial Model"
mean sd
size for nbinomial_1 zero-inflated observations 0.56 0.10
zero-probability parameter for zero-inflated nbinomial_1 0.06 0.06
Scores
Score Poisson Nbin
1 DIC 1111.36 795.44
2 Mlik -579.97 -424.72
3 WAIC 1207.12 798.86
\[ \begin{eqnarray} \text{Model 1: }\ & P(Y=0) &= \pi + (1-\pi)\ \text{Poisson}(\lambda)\\ \text{Model 2: }\ & P(Y=0)&= \pi + (1-\pi)\ \text{NegBinom}(n,p)\\ \end{eqnarray} \] We check the distribution for one person, non-fisher and no children
df_pred = data.frame(nofish = 0, persons = 1, child = 0)
prob_pois = predict(fit_pois, df_pred,
formula = ~ {
pi0 <- zero_probability_parameter_for_zero_inflated_poisson_1
lambda = exp(Intercept + fishing + persons + child)
yy = 0:20
list(
prob0 = pi0 * c(1,rep(0,20)) + (1-pi0) * dpois(yy, lambda = lambda))
} )
prob_nbin = predict(fit_nbin, df_pred,
formula = ~ {
pi0 <- zero_probability_parameter_for_zero_inflated_nbinomial_1
size = exp(size_for_nbinomial_1_zero_inflated_observations)
mu = exp(Intercept + fishing + persons + child)
p = size/(size + mu)
yy = 0:20
list(prob0 = pi0 * c(1,rep(0,20)) + (1-pi0) * dnbinom(yy,size = size, mu = mu))
})Note: To get the names to input in predict use the bru_names() function:
inlabruThere are two ways to implement hurdle models in inlabru
Strategy 1 Use a modified version of the same distributions defined for the zero-inflated models
Strategy 2 Use a “two-likelihood” trick
Available models (Type 0)
Poisson (zeroinflatedpoisson0)
Binomial (zeroinflatedbinomial0)
Negative Binomial (zeroinflatednbinomial0)
BetaBinomial (zeroinflatedbinomial0)
Advantages
It works exactly as the zero-inflated models
No need for extra coding
Disadvantages
Can only be used for the models that are already implemented
The probability of zero (\(\pi\)) is fixed and cannot have covariates
The Model
\[ \begin{aligned} \text{Model 1 }: y_i|\eta_i & \sim \pi1_{(y=0)} + (1-\pi)\text{TruncPois}(\lambda(\eta_i))\\ \text{Model 2 }: y_i|\eta_i & \sim \pi1_{(y=0)} + (1-\pi)\text{TruncNegBin}(n(\eta_i),p(\eta_i))\\ \text{Linear Predictor }: \eta_i & = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} \end{aligned} \]
The code
# define model components
cmp = ~ Intercept(1) + fishing(nofish, model = "linear") +
persons(persons, model = "linear") +
child(child, model = "linear")
# define model predictor
formula = count ~ .
# build the observation model
# Poisson model
lik_pois = bru_obs(formula,
family = "zeroinflatedpoisson0",
data = zinb)
# Negative binomial model
lik_nbin = bru_obs(formula,
family = "zeroinflatednbinomial0",
data = zinb)
# fit the model
fit_pois0 = bru(cmp, lik_pois)
fit_nbin0 = bru(cmp, lik_nbin)Advantages
More flexible model
Can model also the probability of zero \(\pi\)
Can also be used for continuous distributions with mass at 0
Disadvantage
The Hurdle Model handles this by splitting the data-generating process into two parts:
Zero-hurdle component: This part is a binary model that estimates the probability of a zero count versus a non-zero count.
Positive component: This part models the distribution of the positive (non-zero) values.
This idea can be used in inlabru by using two likelihood
A binomial likelihood to model the 0/1 process
A continuous positive density or a truncated count distribution for the positive part
Daily precipitation in Parana state
\[ z_{st} = \left\{ \begin{eqnarray} 0 & \text{ if } &\text{ no rain on day } t \text{ at station } s\\ 1 & \text{ if } &\text{ rain on day } t \text{ at station } s \end{eqnarray} \right. \]
\[ y_{st} = \left\{ \begin{eqnarray} \text{NA} \text{ if } && \text{ no rain on day } t \text{ at station } s\\ \text{rain amout} \text{ if } && \text{ rain on day } t \text{ at station } s \end{eqnarray} \right. \]
\[ \begin{eqnarray} z_{st} & \sim \text{Binom}(p_{st}, n_{st} = 1);\\ (y_{st}|z_{st}=1) &\sim\text{Gamma}(a,b) \\ E(y_{st}) &= \mu_{st} = a/b \end{eqnarray} \]
\[ \text{Where}: \qquad \begin{eqnarray} \eta_t^1 & = \text{logit}(p_t) & = \beta_0^B+u^B(s) \\ \eta_t^2 & = \log(\mu_t)& = \beta_0^G+u^G(s) \end{eqnarray} \]
Implementation
cmp = ~ -1 + Intercept_z(1) + Intercept_y(1) +
space_z(geometry, model = spde) +
space_y(geometry, model = spde) +
local_z(station, model = "iid")
lik_z = bru_obs(formula = z ~ Intercept_z + space_z + local_z,
data = df %>% filter(time==1),
family = "binomial")
lik_y = bru_obs(formula = y ~ Intercept_y + space_y,
data = df %>% filter(time==1),
family = "gamma")
fit = bru(cmp, lik_z, lik_y)Results
We can:
use other components in both linear predictors \(\eta^1\) and \(\eta^2\)
have shared components between the two linear predictor
use other positive-continuous distribution instead of the Gamma (for example log-Normal, Weibull, etc.)
Let’s look again at the fishes examples. We now want to include covariates also in the model for the zero probability:
\[
\begin{eqnarray}
P(Y = 0) = &\pi\\
P(Y=y) = & (1-\pi) \frac{f_Y(y)}{f_Y(0)},& \qquad y =1,2,\dots
\end{eqnarray}
\] We use the same “trick” as before and we now have \[
\begin{eqnarray}
z\sim\text{Binomial}(\pi,n); &\qquad \eta_1 = \text{logit}(\pi) = \beta^1X\\
z\sim\text{Truncated Poisson}(\lambda); &\qquad \eta_2 =\log(\lambda) = \beta^2X\\
\end{eqnarray}
\] In inlabru the truncated Poisson distribution is called nzpoisson
Implementation
cmp = ~ -1 + Intercept_z(1) + Intercept_y(1) +
nofish_z(nofish, model = "linear") + persons_z(persons, model = "linear") +
nofish_y(nofish, model = "linear") + persons_y(persons, model = "linear")
lik_z = bru_obs(formula = z~ Intercept_z + nofish_z + persons_z,
data = zinb,
family = "binomial")
lik_y = bru_obs(formula = y~ Intercept_y + nofish_y + persons_y,
data = zinb,
family = "nzpoisson")
fit = bru(cmp, lik_z, lik_y)Results
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept_z 0.70 0.33 0.04 0.70 1.35 0.70 0
nofish_z 0.23 0.28 -0.33 0.23 0.78 0.23 0
persons_z -0.18 0.12 -0.41 -0.18 0.04 -0.18 0
Intercept_y 0.35 0.15 0.05 0.35 0.65 0.35 0
nofish_y -0.83 0.12 -1.07 -0.83 -0.59 -0.83 0
persons_y 0.53 0.04 0.44 0.53 0.61 0.53 0
In inlabru there is no implementation for the truncated binomial (or negative binomial) distribution
We can “trick” inlabru by using one of the implemented Type 0 models \[
y_i|\eta_i \sim \pi1_{(y=0)} + (1-\pi)\frac{f(y)}{f(0)}
\] with a fixed and small \(\pi\) thus creating a truncated distribution!
Implementation
# Define response for binary and truncated neg. binomial
zinb = zinb %>% mutate(z = ifelse(count==0,1,0),
y = ifelse(count==0,NA, count))
# define model components
cmp = ~ -1 + Intercept_z(1) + Intercept_y(1) +
fishing_z(nofish, model = "linear") + fishing_y(nofish, model = "linear") +
persons_z(persons, model = "linear") + persons_y(persons, model = "linear") +
child_z(child, model = "linear") + child_y(child, model = "linear")